SWSP Project - WESAD Dataset EDA¶

Introduction¶

In this notebook we will explore the extracted features from the WESAD dataset.

%reload_ext pretty_jupyter

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from pretty_jupyter.helpers import matplotlib_fig_to_markdown

pd.set_option('display.max_columns', 100)
pd.set_option('display.max_rows', 100)

Now we will import the data:

data = pd.read_csv('data/m14_merged.csv')
data.head()

data.shape

(1178, 59)

data.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 1178 entries, 0 to 1177
Data columns (total 59 columns):
 #   Column              Non-Null Count  Dtype  
---  ------              --------------  -----  
 0   Unnamed: 0          1178 non-null   int64  
 1   net_acc_mean        1178 non-null   float64
 2   net_acc_std         1178 non-null   float64
 3   net_acc_min         1178 non-null   float64
 4   net_acc_max         1178 non-null   float64
 5   EDA_phasic_mean     1178 non-null   float64
 6   EDA_phasic_std      1178 non-null   float64
 7   EDA_phasic_min      1178 non-null   float64
 8   EDA_phasic_max      1178 non-null   float64
 9   EDA_smna_mean       1178 non-null   float64
 10  EDA_smna_std        1178 non-null   float64
 11  EDA_smna_min        1178 non-null   float64
 12  EDA_smna_max        1178 non-null   float64
 13  EDA_tonic_mean      1178 non-null   float64
 14  EDA_tonic_std       1178 non-null   float64
 15  EDA_tonic_min       1178 non-null   float64
 16  EDA_tonic_max       1178 non-null   float64
 17  BVP_mean            1178 non-null   float64
 18  BVP_std             1178 non-null   float64
 19  BVP_min             1178 non-null   float64
 20  BVP_max             1178 non-null   float64
 21  TEMP_mean           1178 non-null   float64
 22  TEMP_std            1178 non-null   float64
 23  TEMP_min            1178 non-null   float64
 24  TEMP_max            1178 non-null   float64
 25  ACC_x_mean          1178 non-null   float64
 26  ACC_x_std           1178 non-null   float64
 27  ACC_x_min           1178 non-null   float64
 28  ACC_x_max           1178 non-null   float64
 29  ACC_y_mean          1178 non-null   float64
 30  ACC_y_std           1178 non-null   float64
 31  ACC_y_min           1178 non-null   float64
 32  ACC_y_max           1178 non-null   float64
 33  ACC_z_mean          1178 non-null   float64
 34  ACC_z_std           1178 non-null   float64
 35  ACC_z_min           1178 non-null   float64
 36  ACC_z_max           1178 non-null   float64
 37  Resp_mean           1178 non-null   float64
 38  Resp_std            1178 non-null   float64
 39  Resp_min            1178 non-null   float64
 40  Resp_max            1178 non-null   float64
 41  0_mean              1178 non-null   float64
 42  0_std               1178 non-null   float64
 43  0_min               1178 non-null   float64
 44  0_max               1178 non-null   float64
 45  BVP_peak_freq       1178 non-null   float64
 46  TEMP_slope          1178 non-null   float64
 47  subject             1178 non-null   int64  
 48  label               1178 non-null   int64  
 49  age                 1178 non-null   int64  
 50  height              1178 non-null   int64  
 51  weight              1178 non-null   int64  
 52  gender_ female      1178 non-null   int64  
 53  gender_ male        1178 non-null   int64  
 54  coffee_today_YES    1178 non-null   int64  
 55  sport_today_YES     1178 non-null   int64  
 56  smoker_NO           1178 non-null   int64  
 57  smoker_YES          1178 non-null   int64  
 58  feel_ill_today_YES  1178 non-null   int64  
dtypes: float64(46), int64(13)
memory usage: 543.1 KB

We can observe that the all the data is numeric and there are no missing values. We will remove the first column as it is just a clone of the index.

data = data.drop([data.columns[0]], axis=1)
data.head()

Now we will explore the data. We will start by looking at the distribution of the features.

data.describe()

def plot_distribution(data, nbrRows=12, nbrCols=5, figsize=(20, 30)):
    plt.figure(figsize=figsize)
    for c in range(len(data.columns)-12):
        plt.subplot(nbrRows, nbrCols, c+1)
        plt.hist(data.iloc[:, c], bins=50)
        plt.title(data.columns[c])
    plt.tight_layout()
    plt.show()

plot_distribution(data)

We can observe that the net_acc, EDA_phasic and EDA_smna features are very skewed to the left and there may be some outliers. Furthermore observing the acceleraton features we can observe that the readings from the three axis are very similar and it would be a good idea to explore it further. In order to get a better idea how the features correlate with each other we will plot a correlation matrix. Before we plot the matrix

# plot a correlation matrix
corr = data.corr()
mask = np.triu(np.ones_like(data.corr()))
plt.figure(figsize=(40, 40))
sns.heatmap(corr, annot=True, fmt='.2f', cmap='coolwarm', mask=mask)
plt.show()

We can observe that some parts are more saturated than others, these spots indicate features with higher correlation to each other, thus we will explore them further. We will generally ignore high correlation between standard variation, minimum and maximum values as they are expected to be highly correlated. Keeping these factors in mind we can conclude that:

• All the EDA features are highly correlated to each other and the other variables, except for the EDA_smna_min
• The net_acc features are highly correlated to the other variables
• The ACC features are highly correlated to the other variables and we can confirm, our previous observation, that they are also highly correlated to each other.
• BVP and RESP have a very high correlation to other features, with the exception of the mean of those features which have almost no correlation at all.
• TEMP min, max and mean features have a low correlation to the other variables and are incredibly similar to each other, thus they can be combined into one feature. On the other hand the standard deviation and the slope of the temperature have a very low coorelation score with the other variables, thus they can be removed.

For the categorical features the only useful correlations are between the age, height and weight features and the physiological and motion features, thus we can safely remove the rest. Before we do that, it must be noted that there is an obvious correlation between subject and TEMP features, gender and acceleartion, sport and acceleration and of course between gender and height and weight.

Now we will remove the features that won't be of use to us or overlap with other features. The columns in question are: EDA_smna_min, BVP_mean, TEMP_std, TEMP_slope, RESP_mean, smoker_YES, smoker_NO, feel_ill_today_YES, sport_today_YES, coffee_today_YES, gender_male and gender_female. Furthermore we will leave only the TEMP_meanfeature and remove TEMP_min and TEMP_max, in adddition to merging all the acceleration features into one.

columns_to_remove = ['EDA_smna_min', 'BVP_mean', 'TEMP_std', 'TEMP_slope', 'Resp_mean', 'smoker_YES', 'smoker_NO', 'feel_ill_today_YES', 'sport_today_YES', 'coffee_today_YES', 'gender_ male', 'gender_ female', 'TEMP_min', 'TEMP_max']
data_m1 = data.drop(columns_to_remove, axis=1)
data_m1.head()

In order to decide what to do with the acceleration values will again make a correlation matrix, but this time only for the acceleration features.

corr = data_m1.loc[:, 'ACC_x_mean':'ACC_z_max'].corr()
mask = np.triu(np.ones_like(corr))

plt.figure(figsize=(10, 10))
sns.heatmap(corr, annot=True, fmt='.2f', cmap='coolwarm', mask=mask)
plt.show()

As we can see th columns are highly correlated to each other, thus we will merge the ACC_x, ACC_y and ACC_z features into ACC_mean_mean, ACC_mean_std, ACC_mean_min, ACC_mean_max.

# combine the ACC features for the three axis into one feature
data_m1['ACC_mean_mean'] = data_m1[['ACC_x_mean', 'ACC_y_mean', 'ACC_z_mean']].mean(axis=1)
data_m1['ACC_mean_std'] = data_m1[['ACC_x_std', 'ACC_y_std', 'ACC_z_std']].mean(axis=1)
data_m1['ACC_mean_min'] = data_m1[['ACC_x_min', 'ACC_y_min', 'ACC_z_min']].mean(axis=1)
data_m1['ACC_mean_max'] = data_m1[['ACC_x_max', 'ACC_y_max', 'ACC_z_max']].mean(axis=1)
data_m2 = data_m1.drop(['ACC_x_mean', 'ACC_y_mean', 'ACC_z_mean', 'ACC_x_std', 'ACC_y_std', 'ACC_z_std', 'ACC_x_min', 'ACC_y_min', 'ACC_z_min', 'ACC_x_max', 'ACC_y_max', 'ACC_z_max'], axis=1)
data_m2.head()

!!!!!!!!!!!!!!!!!!!! REMOVE OR REVISIT !!!!!!!!!!!!!!!!!!!!!!!
We will also temporarily remove the 0_mean, 0_std, 0_min and EDA_max features as they do not make any sense.

data_m2 = data_m2.drop(['0_mean', '0_std', '0_min', '0_max'], axis=1)

Now we will plot the correlation of the features to the target variable.

# one hot encode the label column
data_m2_ohe = pd.get_dummies(data_m2, columns=['label'], prefix='', prefix_sep='')

corr = data_m2_ohe.corr()
corr_data = pd.DataFrame({'amusement': corr['0'], 'baseline': corr['1'], 'stress': corr['2']}) 
corr_data = corr_data.drop(['0', '1', '2'], axis=0).sort_values(by='baseline', ascending=False)

fig, ax = plt.subplots(figsize=(10,10))
sns.heatmap(corr_data, vmin=-1, vmax=1, cmap="seismic", annot=True, ax=ax)
plt.show()

From the graph above we can conclude that the feature that have the highest correlation to the target vriable are EDA_tonic, Resp and net_acc.

data_m2.label.unique()

array([1, 2, 0], dtype=int64)

We can see that there are three unique values in the label column, according to the WESAD data preparation repository these values represent the following stress levels:

• 0 - Amusement
• 1 - Baseline
• 2 - Stress

Now we will separate the data by each label and observe how the features are affected by different stress levels.

amusement_data = data_m2[data_m2['label'] == 0].reset_index(drop=True)
baseline_data = data_m2[data_m2['label'] == 1].reset_index(drop=True)
stress_data = data_m2[data_m2['label'] == 2].reset_index(drop=True)

def plot_data_per_label(data, feature, labels=['amusement', 'baseline', 'stress']):	
    plt.figure(figsize=(20, 6))
    plt.title(f"{feature} variance per label")
    plt.xlabel("N records (30 second window)")
    plt.ylabel("{feature} measurement")

    for d in data:
        plt.plot(d[feature])

    plt.legend(labels, loc='upper right')
    plt.show()

The features that we will plot will be the same as the ones that we used to plot the correlation to the target variable, except the subject, age, height and weight as they had a very low correlation to the label, thus we won't observe any meaningful difference in the plots.

features = data_m2.columns.drop(['label', 'subject', 'age', 'height', 'weight'])
features

Index(['net_acc_mean', 'net_acc_std', 'net_acc_min', 'net_acc_max',
       'EDA_phasic_mean', 'EDA_phasic_std', 'EDA_phasic_min', 'EDA_phasic_max',
       'EDA_smna_mean', 'EDA_smna_std', 'EDA_smna_max', 'EDA_tonic_mean',
       'EDA_tonic_std', 'EDA_tonic_min', 'EDA_tonic_max', 'BVP_std', 'BVP_min',
       'BVP_max', 'TEMP_mean', 'Resp_std', 'Resp_min', 'Resp_max',
       'BVP_peak_freq', 'ACC_mean_mean', 'ACC_mean_std', 'ACC_mean_min',
       'ACC_mean_max'],
      dtype='object')

We can see that the categorical values have been successfully removed, thus we will proceed to plot the efect of the stress level on the features.

Feature variance per label¶

.tabset|O_O|.tabset-pills|O_O|.tabset-fade

from scipy.interpolate import make_interp_spline

# required inputs: data - a list of dataframes that correspond to the data for each label, feature - the feature to plot, labels - the labels for the data
def plot_smooth_data_per_label(data, feature, labels=['amusement', 'baseline', 'stress'], num=100, b_spline_degree=3):
    plt.figure(figsize=(20, 6))
    plt.title(f"{feature} variance per label")
    plt.xlabel("N records (30 second window)")
    plt.ylabel(f"{feature} measurement")

    # loop through the data for each srtress level, smooth it and plot it
    for d in data:
        xnew = np.linspace(d[feature].index.min(), d[feature].index.max(), num)
        spl = make_interp_spline(d[feature].index, d[feature], k=b_spline_degree)
        power_smooth = spl(xnew)
        plt.plot(xnew, power_smooth)

    plt.legend(labels, loc='upper right')
    #plt.show()
    return plt.gcf()

stress_level_data = [amusement_data, baseline_data, stress_data]
#for f in features:
#    plot_smooth_data_per_label(stress_level_data, f, num=50, b_spline_degree=2)

from IPython.display import display, Markdown
for f in features:
    display(Markdown(f"## {f} feature"))
    plot = plot_smooth_data_per_label(stress_level_data, f, num=50, b_spline_degree=2)
    display(Markdown(matplotlib_fig_to_markdown(plot)))

net_acc_mean feature¶

net_acc_std feature¶

net_acc_min feature¶

net_acc_max feature¶

EDA_phasic_mean feature¶

EDA_phasic_std feature¶

EDA_phasic_min feature¶

EDA_phasic_max feature¶

EDA_smna_mean feature¶

EDA_smna_std feature¶

EDA_smna_max feature¶

EDA_tonic_mean feature¶

EDA_tonic_std feature¶

EDA_tonic_min feature¶

EDA_tonic_max feature¶

BVP_std feature¶

BVP_min feature¶

BVP_max feature¶

TEMP_mean feature¶

Resp_std feature¶

Resp_min feature¶

Resp_max feature¶

BVP_peak_freq feature¶

ACC_mean_mean feature¶

ACC_mean_std feature¶

ACC_mean_min feature¶

ACC_mean_max feature¶

Findings from feature variance¶

From the plots above we can clearly see that in general the net_acc_std, RESP_max, RESP_std and all the EDA features are higher when a person is stressed. The ACC, TEMP, Resp_min and BVP features do not show any sustantial difference between the stress levels. Additionally, there is not any absorvabale difference between baseline and amusement stress levels that is correlated to the features.

Now we will use the SelectKBest method to select the best features. We will use the f_classif method as we are dealing with numerical input and categorical output.

from sklearn.feature_selection import SelectKBest, f_classif
empE4_features = features.to_list()
empE4_features.remove('Resp_std')
empE4_features.remove('Resp_min')
empE4_features.remove('Resp_max')
data_m3_X = data_m2[empE4_features]
data_m3_y = data_m2['label']
print(f"Old data shape: {data_m3_X.shape}")	

data_m3_X_new = SelectKBest(f_classif, k=5).fit_transform(data_m3_X, data_m3_y)
print(f"New data shape: {data_m3_X_new.shape}")

selected_features = data_m3_X.columns[SelectKBest(f_classif, k=5).fit(data_m3_X, data_m3_y).get_support()]
selected_features

Old data shape: (1178, 24)
New data shape: (1178, 5)

Index(['net_acc_std', 'net_acc_max', 'EDA_tonic_mean', 'EDA_tonic_min',
       'EDA_tonic_max'],
      dtype='object')

We can observe that the selected features are 'net_acc_std', 'net_acc_max', 'EDA_tonic_mean', 'EDA_tonic_min', 'EDA_tonic_max', 'Resp_std', 'Resp_max'. Now we will plot each feature against the target variable per subject to observe if there is a clearer diferentiation between the stress levels in the selected features .

Data analysis by subject¶

.tabset|O_O|.tabset-pills|O_O|.tabset-fade

def plot_smooth_data_per_label_per_subject(data, feature, subjects, labels=['amusement', 'baseline', 'stress'], num=100, b_spline_degree=3):
    plt.figure(figsize=(20, 20))

    for s in subjects:
        s_index = np.where(subjects == s)[0][0]
        plt.subplot(5, 3, s_index + 1)
        plt.title(f"{feature} variance per label for subject {s}")
        plt.xlabel("N records (30 second window)")
        plt.ylabel(f"{feature} measurement")
        
        data1 = [data[0][data[0]['subject'] == s].reset_index(drop=True), 
            data[1][data[1]['subject'] == s].reset_index(drop=True), 
            data[2][data[2]['subject'] == s].reset_index(drop=True)]

        # loop through the data for each stress level, smooth it and plot it
        for d in data1:
            xnew = np.linspace(d[feature].index.min(), d[feature].index.max(), num)
            spl = make_interp_spline(d[feature].index, d[feature], k=b_spline_degree)
            power_smooth = spl(xnew)
            plt.plot(xnew, power_smooth)
        plt.legend(labels, loc='upper right')

    plt.subplots_adjust(left=0.1,
                    bottom=0.1,
                    right=1,
                    top=0.9,
                    wspace=0.3,
                    hspace=0.4)
    return plt.gcf()

data_m3 = data_m2[[*selected_features.tolist(), 'label', 'subject']]
amusement_data = data_m3[data_m3['label'] == 0].reset_index(drop=True)
baseline_data = data_m3[data_m3['label'] == 1].reset_index(drop=True)
stress_data = data_m3[data_m3['label'] == 2].reset_index(drop=True)

#for f in selected_features.tolist():
#    plot_smooth_data_per_label_per_subject([amusement_data, baseline_data, stress_data], feature=f, subjects=data_m2['subject'].unique())

for f in selected_features.tolist():
    display(Markdown(f"## {f} variance per label for all subjects"))
    plot = plot_smooth_data_per_label_per_subject([amusement_data, baseline_data, stress_data], feature=f, subjects=data_m2['subject'].unique())
    display(Markdown(matplotlib_fig_to_markdown(plot)))

net_acc_std variance per label for all subjects¶

net_acc_max variance per label for all subjects¶

EDA_tonic_mean variance per label for all subjects¶

EDA_tonic_min variance per label for all subjects¶

EDA_tonic_max variance per label for all subjects¶

Conclusion¶

From observing the graphs above we can conclude that for the EDA_tonic features almost all the subjects, except for subject 17, have a noticeably higher and fluctoating values when they are stressed. In the case of the net_acc_std feature we can see that there are way more fluctiotaions in the readings of stress patients and the values are generally higher, the difference in the values of stressed and baseline periods can be observed more clearly in the net_acc_max feature graphs as there aren't any big oscilations. We can conclude that even a basic Machine Learning model, or a thresholding algorithm will be able to classify the stress levels with a high accuracy, because of the clear difference between the features of the different stress levels.